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Abstract 

<N : 

We review some aspects of current knowledge regarding the decay of metastable 
f-*>) . phases in many-particle systems. In particular we emphasize recent theoretical and 

computational developments and numerical results regarding homogeneous nucle- 
ation and growth in kinetic Ising and lattice-gas models. An introductory dis- 
cussion of the droplet theory of homogeneous nucleation is followed by a discus- 
sion of Monte Carlo and transfer-matrix methods commonly used for numerical 
study of metastable decay, including some new algorithms. Next we discuss specific 
classes of systems. These include a brief discussion of recent progress for fluids, 
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and more exhaustive considerations of ferromagnetic Ising models (i.e., attractive 
lattice-gas models) with weak long-range interactions and with short-range inter- 
actions. Whereas weak-long-range-force (WLRF) models have infinitely long-lived 
metastable phases in the infinite-range limit, metastable phases in short-range-force 
(SRF) models eventually decay, albeit extremely slowly. Recent results on the finite- 
size scaling of metastable lifetimes in SRF models are reviewed, and it is pointed 
out that such effects may be experimentally observable. 



1 Introduction 

Metastable phases are common in nature, and some of them have such extremely long 
lifetimes that they are practically indistinguishable from equilibrium phases. For example, 
even though the stable phase of carbon at room temperature and atmospheric pressure 
is graphite, the claim that diamonds are forever is not likely to be challenged on physical 
grounds. Other examples of metastable phases with average lifetimes that are measured 
in milliseconds to days, rather than in millions or billions of years, are supercooled or 
supersaturated fluids [1-32], ferroelectrics [33-38], the small single-domain ferromagnetic 
particles important in paleomagnetism and high-density recording media [39-49], and 
vortex states in superconductors f5"0fl . Further possible examples are the supercooled 
quark/gluon plasma associated with the QCD confinement transition |51], |52|] and the 
"false vacuum" associated with the electroweak transition [53-57], both of which may 
have played important roles in the early development of the universe. 

To clarify what we understand by a metastable phase, it is hard to improve on the 



empirical descriptions due to Penrose and Lebowitz [|58], ^] and Sewell |K) . 
(i) The free energy of the system is not fully minimized. 

(ii) Only one thermodynamic phase is present, and for sufficiently small and slow pertur- 
bations the usual laws of reversible thermodynamics apply. 

(Hi) A system starting in the metastable phase is likely to take a long time to escape. 
(iv) Escape is irreversible: once out of the metastable phase, the system is extremely 
unlikely to return. 

We concentrate on systems for which the order parameter is a nonconserved scalar 
(Model A in the Hohenberg-Halperin scheme of dynamic universality classes ||61|| ) so that 
the metastability is imposed by an applied external field. (The closely related phenomenon 
of hysteresis imposed by an oscillating field is discussed by Acharyya and Chakrabarti 
elsewhere in this volume [j£>2|| .) Even so, many of the phenomena we discuss can be 
generalized to systems with a multidimensional and/or conserved order parameter . 

The wide variety of contexts in which metastability has been studied, often indepen- 
dently, makes it difficult to present a discussion that is equally accessible to all potentially 
interested readers. In this review we have chosen to use the language of Ising models. 
The most important reason for this choice is our belief that the present theoretical un- 
derstanding of metastable phases and their modes of decay is most highly developed in 
applications to the relatively simple kinetic Ising models and their equivalent formulations 
as lattice-gas models. A secondary reason is the high symmetry of the Ising formulation. 
Throughout this paper we therefore use Ising ferromagnets as generic examples. Below we 
define these models and give the standard transformations that will enable the reader to 
convert our results to the lattice-gas language appropriate to, e.g., fluids, binary mixtures, 
and adsorbate systems. 

An Ising ferromagnet is defined by a Hamiltonian (i.e., an energy functional), 

where Sj = ±l ("up" or "down") is a binary variable, or "spin," at site i, H is the applied 
field, and the sums run over all M sites on a <i-dimensional lattice. The interaction energies 
Jij are positive and symmetric under interchange of the site indices i and j, and without 



loss of generality we set J^i—O. By requiring that J2j -h,j = zJ is independent of i we 
ensure the spatial homogeneity of the energy functional. We also note that TC is invariant 
under the transformation {si —>■ —S{, H —* —H}. The order parameter conjugate to H 
is the "magnetization" or "polarization" , 
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■A/^E*. (2) 



As is well known |B3 , this model is equivalent to a two-state lattice-gas model with local 



concentration variables Cj=0 or 1 ("empty" or "occupied"), for which Eq. (|1|) takes the 
form 

ft = -g E ®i,j c * c i ~ t* E c * + y (z 1 - 2^°) - ^ LG + y (^ ~ 2^°) ' ( 3 ) 

where the quantities appearing in the two equivalent formulations of the Hamiltonian, 
Eq. ([1]) and Eq. (|3|), are linked by the transformations 

a = (s t +i)/2 

$. . = 4/. • 
// = 2#+/i . (4) 

Here <J>jj are attractive lattice-gas interaction energies and // is the chemical potential, 
whose value at coexistence (i.e., for H=0) is [Iq=—1zJ. The chemical potential is related 
to the (osmotic) pressure p as //— // = ^B^l n (p/Po); where /cb^ is Boltzmann's constant 
times the absolute temperature, andpo is the pressure at coexistence. The order parameter 
conjugate to /i is the density, 

P = AT 1 $> = (™ + l)/2. (5) 

i 

The energy functional TC is not a quantum-mechanical Hamiltonian and does not im- 
pose a unique dynamic on the system. When one uses an Ising or lattice-gas model to 
describe the kinetics of a physical system, one must therefore define a specific stochastic 
dynamic to simulate the interactions with the system's surroundings. Usually it is phys- 
ically most realistic to choose a dynamic which is local in the sense that it only allows 
transitions involving a single site or a pair of nearest-neighbor sites, and in this review 



we mostly limit our attention to such dynamics. Simple examples are the Metropolis pi | 
and Glauber |)5j dynamics, which are among the algorithms we consider in the context of 
Monte Carlo simulations in Sec. |3.1| . (Under somewhat restrictive conditions, the Glauber 
dynamic has been obtained from the quantum mechanical equations of motion for a sys- 
tem of distinguishable spin- 1/2 particles weakly coupled to an infinitely large quantum 
reservoir J66| .) 

The Ising (lattice-gas) model below its critical temperature and in the absence of an 
applied field (with /x=/i ) has two coexisting, ordered phases of equal free energy: one 
with positive magnetization (high density) and one with negative magnetization (low 
density). This degeneracy is lifted by applying a nonzero field (changing \x away from 
Ho): the phase whose magnetization is parallel to the field (for which p(/j,)—p(no) has the 



same sign as /x— /xo) becomes the unique equilibrium phase, and the one with the opposite 
magnetization ("wrong" density) becomes metastable. 

The system can be prepared in the metastable phase by equilibrating it in a nonzero 
field that is then instantaneously reversed. Although the system is no longer in equilibrium 
immediately after the field reversal, it is nevertheless stable against small fluctuations, 
and its thermodynamic properties are similar to those it would have in the equilibrium 
phase. This is because configurations obtained by flipping small clusters of neighboring 
spins cost more free energy by introducing an interface between previously parallel spins 
than they gain by lowering the field energy. In order for this metastable phase to decay, 
it is therefore necessary that a cluster is created that is sufficiently large for the free 
energy gained by aligning more spins with the field to just outweigh the cost of breaking 
the necessary extra bonds. Such a fluctuation corresponds to a local maximum in the 
free-energy landscape and is usually called a "critical nucleus" or a "critical droplet". 
Once randomly created through a thermal fluctuation in the metastable phase, it is likely 
to grow further. During this growth period the now "supercritical" droplet incorporates 
spins from the metastable phase into a growing, and soon macroscopic, domain of the 
equilibrium phase. The timescale for the creation of a critical droplet is in general much 
longer than that characteristic of the subsequent growth. 

Real systems are of course more complicated than the simple picture discussed above. 
Nevertheless, it covers the essential physics of nucleation in a variety of situations. Some 
of the additional complications that arise in the study of metastability in fluids with 
continuous degrees of freedom are briefly discussed in Sec. £|. 

The nucleation mechanism described above, which is usually known as homogeneous 
or thermally activated nucleation, will be the focus of our attention. It is the dominant 
mechanism in systems that are free of defects, or in which the defect concentration is low, 
or in which the defects are much smaller than the critical droplet size. Whenever this 
is not the case, heterogeneous nucleation on defects or interfaces may be the dominant 



process for producing droplets of the equilibrium phase |HJ. However, the subsequent 
growth process depends little on the nucleation mechanism |33], |68 . 

The organization of the remainder of this review is as follows. Section ^| contains 
a brief review of classical nucleation theory and some of the general results obtained 
in the "post-classical" period after the 1940's. In Sec. |^ we review numerical Monte 
Carlo (in Sec. |3.1|) and transfer-matrix (in Sec. |Q|) methods to study nucleation and 
metastable decay, primarily in kinetic Ising models. Following these general sections we 
consider specific classes of systems: fluids in Sec. |], Ising models with weak long-range 
forces (WLRF models) and their relation to mean-field approximations in Sec. |^, and 
Ising models with short-range forces (SRF models) in Sec. |5|. The subsections in Sec. |5] 



are Sec. 5.1, which contains specific nucleation rates for the ramified critical droplets 



characteristic of WLRF models, and Sec. 5.2, which contains numerical results. The 



subsections in Sec. ^| are Sec. |6.1| , which contains specific nucleation rates for the compact 
critical droplets characteristic of SRF models, Sec. |6.2| , in which we discuss the interplay 



between nucleation and the growth of supercritical droplets, Sec. |673|, in which we consider 
finite-size effects, and Sec. |6^ , in which we consider effects of the discrete lattice at low 
temperatures. Finally, in Sec. |7|, we give a concluding summary and discussion. 



2 Droplet Theory of Homogeneous Nucleation 

The first scientific description of metastability may well be Fahrenheit's experiments with 
supercooled water, published in 1724 IT], an d further observations were reported during 



the eighteenth and early nineteenth centuries by several workers ||19| . The basis for a the- 
oretical understanding was laid 150 years after Fahrenheit's paper with van der Waals' 
and Maxwell's || early mean-field approaches and Gibbs' realization that the reversible 
work of formation of a droplet of the stable phase in a metastable background is a "mea- 
sure of the stability" of the metastable phase [4-6]. Although Gibbs' work clearly states 
the basic equations of what is today known as "Classical Nucleation Theory" (CNT), 
further development did not occur until the period from the 1920's through the 1940's. 
Gibbs' essentially thermodynamic (i.e., static) approach was then complemented by ki- 
netic considerations by Volmer and Weber J7|, Farkas ||, Becker and Doring ||, and 



Zeldovich [|T0|, pT| , whereas Bijl |12|| , Frenkel fl3|l , and Band |I4| focused on calculating 
a constrained, metastable partition function through a cluster-expansion procedure in- 
spired by that of Mayer p9l. The basic result for the nucleation rate can be written in a 



Van't Hoff-Arrhenius |70|, [7TJ form, 

T = Ve~^ , (6) 

where V is a non-exponential prefactor, F c is the free-energy cost of a critical droplet 
(Gibbs' "measure of stability"), and (3=1/ k^T is the inverse temperature. (A compact 
historical sketch is given by Dunning |l9|. For more recent general reviews of metastability 
and the kinetics of first-order phase transitions see, e.g., Refs. [31, 72-74].) 

The "post-classical" developments in nucleation theory are mainly efforts to evaluate 
F c and/or the prefactor in Eq. ((^) for specific systems. This observation in itself provides 
an important insight: the lifetime of a metastable state depends crucially on the structure 
and dynamics of the excitations through which it decays ]75| . 



The relevance of the classification scheme for dynamic critical phenomena due to 
Hohenberg and Halperin |61| to metastable decay is not clear. In addition to those 



aspects of the physics that influence both the static and the dynamic universality, such 
as the interaction range and the symmetries and dimensions of both the system itself and 
the order parameter, the conservation laws that influence the dynamic scaling should also 
be important for the nucleation rate. 

The "metastable thermodynamics" that can be observed during the period before the 
decay takes place is almost indistinguishable from true equilibrium thermodynamics. This 
has inspired efforts to treat metastability through suitable generalizations of equilibrium 
statistical mechanics. Notable among these are three papers by Langer [76-78], in which 
he shows that for a wide class of models, whose dynamics can be described by a Fokker- 
Planck equation, the homogeneous nucleation rate (per unit time and volume) may be 
written as 

r = ^|Im/|, (7) 

where Im/ is the imaginary part of a "metastable" free-energy density /, which is ob- 
tained by analytically continuing the equilibrium free-energy density / into the metastable 
phase. The "kinetic prefactor" k contains all dependence on the specific dynamic. This 



important work brings together aspects of CNT with results on droplet theory and ana- 
lytic continuation from Andreev |79[ and Fisher [RDI and on thermally activated processes 



from Landauer and Swanson 811 and Kramers |8 



Despite its apparent simplicity and its similarity to the corresponding relation in quan- 
tum mechanics (obtained by replacing / by the energy and fin/it by 2/h) no general proof 
of Eq. (^) exists, and its domain of validity remains unclear. Analytic verification for spe- 
cific models has been reported, e.g., by Newman and Schulman [f^| for a Curie- Weiss 
ferromagnet, by Roepstorff and Schulman |53J for an urn model, by Gaveau and Schul- 
man |85[ for a class of models including the droplet model of condensation ||TT], [79|, |8*0|1 , 



and by Penrose [§§] for the droplet model with Becker-Doring dynamics. A different 
approach is to assume the validity of Eq. (|7]) and calculate |Im/| analytically or numer- 
ically. Again, analytic calculation requires a specific model for the fluctuations included 
in the calculation of the analytic continuation /. Such field-theoretical calculations were 
done by Coleman and Callan p5| , 53 for the "false vacuum" in quantum-field theory, 



by Biittiker and Landauer [37], |88| for a one- dimensional (ID) overdamped sine-Gordon 



chain, by McCraw |g9] for a ID Kac model with algebraically decaying interactions, by 
Giinther, Nicole, and Wallace |Xj, who generalized Langer's field-theoretical calculation 
to arbitrary spatial dimension (see Sec. |6 . 1|) , by Zwerger [91] for a 4 field model, by 

31 for 



Cottingham et al. |32j for a model of bubble formation in fluids, by Braun 
a model of switching in single-domain ferromagnetic particles, and by Klein and Unger 
and Gorman et al. |94], |9^ 



93 



who used 3 field theories to study WLRF systems 
near the classical spinodal (see Sec. ^). 

As may be understood from the heuristic discussion of the decay of a metastable phase 
presented in Sec. [I], the metastable lifetime depends not only on the rate of nucleation 
of critical droplets, but equally importantly on the subsequent rate of growth of the 
supercritical droplets and on the interplay between these processes of nucleation and 
growth. This was realized by Kolmogorov |96|| , Johnson and Mehl [97j, and Avrami |98| 



(KJMA) at about the same time as CNT was being developed, and a few years later by 
Evans 1681. 



3 Computational Methods 

In contrast to analytical calculations of F or |Im/|, nonperturbative numerical methods 
do not require a priori knowledge about the critical excitations. In this section we briefly 
discuss two classes of such methods: some varieties of the well-known Monte Carlo (MC) 
method for numerical simulation, and a recently introduced extension of the equilibrium 
transfer-matrix (TM) method to also encompass metastable phases. 



3.1 Monte Carlo Methods 

Simulations using standard Metropolis or heat-bath dynamics with spatially local updates 
(see, e.g., Ref. |99|) have remained the computational methods of choice for studying 
metastable decay in systems with nonconserved order parameter. The spatial locality of 
the algorithm preserves the free-energy barriers that dominate the nucleation rate, leading 
to relatively faithful representations of metastable dynamics in physical systems. (This 



is not generally true for cluster algorithms such as the one due to Swendsen and Wang 
[ |100|| , which are used to accelerate equilibrium simulations near criticality. However, see 
further discussion below.) 

In equilibrium studies, the choice between the Metropolis |64| and heat-bath algo- 



rithms (the latter also known as the Gibbs sampler or, when applied to the Ising model, 
the Glauber ^5j dynamic) is largely a matter of convenience, and it is often considered 
so trivial that it is not even stated explicitly. This habit sometimes has carried over to 
dynamical studies where, in our opinion, attention should be given to finding the dynamic 
best representing the physical system to be simulated. For completeness we give below 
the probabilities for an allowed transition from state x to x' for these two algorithms 
[101-103]. 

In the Metropolis algorithm a candidate state is selected according to a proposal dis- 
tribution, which vanishes for forbidden transitions. (The decision about which transitions 
should be allowed is part of the full definition of the algorithm. For an Ising system 
one could for instance allow single-spin flips, nonlocal cluster flips as in the Swendsen- 



Wang ||100| and related [104-106] algorithms, or nearest-neighbor spin exchanges as in the 
conserved-order-parameter Kawasaki dynamic [|107|| .) The proposed transition is accepted 
with probability one if it leads to a reduction in energy, whereas an increase in energy is 
accepted with probability given by a Boltzmann factor, so that the acceptance probability 
may be written as 

W M (x^x') = min{l,exp[-P(E(x')-E(x))]} . (8) 

In contrast, the heat-bath algorithms accept any state x' to which a transition from x is 
allowed, with the equilibrium probability over the set of all accessible states x": 

W R (x -> x') = exp [-(3E{x')\ I J2 ex P [-&&&')] • (9) 

/ {x" accessible from x} 

Both the Metropolis and the heat-bath algorithm satisfy detailed balance, and with 
ergodicity they eventually converge to thermodynamic equilibrium. However, the detailed 
Markov processes generated are not identical. Beside the choice of the Metropolis or heat- 
bath transition probability, more subtle differences between algorithms may also influence 
the results of simulations, and even though data obtained by different local algorithms 
can often be connected by simply rescaling the time, the rescaling factor may depend 
nontrivially on temperature and distance from the coexistence curve. For example, in a 
recent study of metastable decay in the two-dimensional Ising model with the Metropolis 
dynamic ||108|| , Rikvold et al. found that the manner in which the candidate site for 



the next spin update was chosen (sequentially or randomly) affected the observed field 
dependence of the kinetic prefactor in Eq. (|7]). Similarly, in a recent study of magnetization 
relaxation in the three-dimensional Ising model at the critical point, Ito demonstrated that 
the detailed nature of the finite-size effects depends not only on the choice between the 
Metropolis and heat-bath dynamics, but also on whether a sequential or a checker-board 
update scheme was used ||109| . In choosing the MC algorithm for a particular study, 



one should therefore consider whether the quantities of interest are universal in the sense 
that they are the same for all local algorithms, or whether they depend on the dynamical 
details of the algorithm used. 



For low temperatures and close to coexistence, both the local algorithms discussed 
above and the physical dynamic they simulate spend most of the time creating short- 
lived, microscopic excitations in the metastable phase. Because of the smallness of the 
subcritical clusters, this is also true for the modified Swendsen-Wang cluster dynamic 



used in MC simulations by Ray, Tamayo, and Wang ||104j , |105|| and studied theoretically 
by Martinelli, Olivieri, and Scoppola [ |106|| . The brute-force way around this problem is 
to apply more computer power in the form of various kinds of supercomputers or special- 
purpose machines (as, e.g., in Refs. ||108| , |110| , |111| ), but inevitably the slowness of the 
dynamic makes this approach impracticable. Very recently, two novel MC algorithms have 
been introduced, which make simulations deep in the metastable region feasible. Both 
methods utilize absorbing Markov chains ||112| . One of them, introduced by Novotny [113- 
115], generalizes the rejection- free n-fold way algorithm of Bortz, Kalos, and Lebowitz 
"THl and achieves CPU-time savings of many orders of magnitude relative to the local 



algorithms. This speedup is obtained without changing the underlying dynamic in any 
way. The other method, introduced by Lee et al. [ |117|| , combines absorbing Markov chains 
with the Multicanonical method 



11* 



The resulting dynamic depends on the microscopic 
configurations only through their projection onto the order parameter. Nevertheless, most 
qualitative and quantitative features of the dependences of the metastable lifetimes on 
field, temperature, and system size agree with theoretical results and direct simulations for 
local dynamics. By comparison with such more traditional methods, this algorithm may 
help deepen our understanding of universality in metastable decay. Preliminary results 
from both methods are promising (some are presented in Sees. |6^ and |6.4|) , and we hope 
to review further progress in the future. 

It is common in dynamical MC simulations of metastable decay (regardless of the 
particular algorithm employed) to study the relaxation of the order parameter, starting 
from the metastable phase. This approach is closely related to the use of nonequilibrium 
relaxation functions, introduced by Binder [|119|| . It has been used for SRF models in two 
[37, 104, 108, 110, 111, 120-129], three [38, 105, 130-132], and higher g3| dimensions, 
and also for WLRF models p|, |3|, [133, |135| 
Sec. 1. 



Some recent results are presented in 



3.2 The Constrained-Transfer-Matrix Method 



The Constrained- Transfer-Matrix (CTM) method introduced by one of us | 136 [ is partic- 
ularly suited for numerical calculation of the analytically continued free-energy density 
/ in the field-theoretical droplet theory discussed in Sec. |2], without needing a theoret- 
ical model of the critical droplet geometry. The method extends the usual concept of 
the transfer matrix (TM) ||137|| to also include constrained nonequilibrium states. Here 
we give a brief description of the technique. More detailed discussions can be found in 
Refs. H (J|, m |13|, [139 . 

In a standard TM calculation, an NxL lattice is considered in the limit L—kx>, and 
the Hamiltonian (the energy functional) is written as a sum of layer Hamiltonians, 7i = 
Ya=i T^(xi, xi + i), where H depends only on the configurations xi and xi + i of two adjacent 
layers. The TM is a positive matrix, 

T = Y / \x)e- mx ' x ' ) - 



x 



^|a)A c 



{a\ 



(10) 



where the second equality represents a standard eigenvalue expansion. By the Perron- 
Frobenius theorem |137|| , the dominant eigenvalue Aq is positive and nondegenerate, and 



the corresponding eigenvector |0) is the only one with elements that can all be chosen 
positive. From To one can calculate by standard methods the probability densities, cor- 
relation functions, and partition function that fully describe the equilibrium phase ||137|| . 
The CTM method generalizes this well-known equilibrium technique by associating 
with each eigenvalue X a a "constrained" TM T a so that "constrained" joint and marginal 
probability densities are defined in analogy with the equilibrium (a=0) case: 

P a (x u xi +k ) = {a\xi){xi\(\~ 1 T a )^\xi +k ){xi +k \a) 

P a {xi) = (a\xi)(xi\a) . (11) 



It was pointed out by McCraw, Schulman, and Privman ||127|, |14U|, |141|| that the con- 



strained marginal probability densities P a {x), as defined above, can be interpreted as 
actual probability densities over single-layer configurations in a constrained phase. To 
obtain explicitly the constrained joint probability densities P a (xi,xi + k) that define the 
inter-layer correlations in the constrained phase, one must determine the constrained TM, 
T a . It is chosen to commute with T , and in order to ensure convergence of P a (xi, xi + k) 
towards stochastic independence as | A:| — ^oo, each eigenvalue |A|>|A Q | is reweighted to be- 
come A a /A, so that the dominant eigenvalue of T Q is A Q . A "constrained" free-energy 
density is defined by 

In A Q 1 ^ , . / \Xi \T a \Xi+i) \ , . 

f - = —mr + w £, PJx " x,+l)Ln I wf^y j • < 12 > 

where Ln(z) is the principal branch of the complex logarithm. For a=0 this reduces to 
the standard equilibrium result, f — fo — — (NP)" 1 In A . For a>0 the eigenvector |a), 
corresponding to the dominant eigenvalue of T a , is orthogonal to |0) and therefore cannot 
have all positive elements. As a result, T a cannot in general be a positive matrix, and 
the second term in Eq. ( |T2D becomes complex-valued. It has been observed by Giinther 
et al. [ |139|| that this second term can be considered as a complex generalization of the 
Kullback discrimination function (see, e.g., Ref. [|142|| ) for P a (xi,xi+i) with respect to the 



divergent 'probability density' obtained by substituting T for T a in Eq. (|TT|). 

Successful applications of this method to calculate |Im/| have recently been performed, 
both for WLRF systems R94] |95|, |T36J [1431, and for SRF systems |3|, |139|, |H|. Some 
results are shown in Sees. || and |^. 

4 Fluid Systems 

As mentioned in Sec. [I], metastable decay in fluid systems, such as condensation of a super- 
saturated vapor or freezing of a supercooled liquid, presents complications not encountered 
in Ising or lattice-gas systems. In this section we mention some of these difficulties and 
give a few references. 

Most of the early theoretical work on metastability that led to the formulation of CNT 
[4-14] was explicitly concerned with fluids, as were the early mean-field approaches of van 
der Waals and Maxwell |§. Extensive reviews can be found in Refs. [17-19]. 



One of the most obvious differences between fluids and lattice-gas systems is that 
fluids have continuous degrees of freedom associated with droplet translation and rotation. 
Their effects were considered by Lothe, Pound, and collaborators [20, ^IJ, who evaluated 
the prefactor in Eq. (^) and obtained an increase in the predicted condensation rate, 
relative to earlier theories, on the order of 10 17 . (See, however, a recent discussion in 
Ref. " 



25 



The long controversy over the Lothe-Pound result emphasizes the importance 
and difficulty of understanding the pre-exponential factors in the nucleation rate. 

An important conceptual problem arises when one tries to formulate a precise droplet 
definition [25-29]. It was already considered by Gibbs, who introduced a "dividing sur- 
face" between the liquid droplet and the vapor, but it does not seem yet to have reached 
a fully satisfactory solution for fluid systems [|I45||. Even in the much simpler Ising model, 



the proper definition of droplets has been realized only relatively recently in terms of the 
"Fortuin-Kastelyn-Coniglio-Klein-Swendsen-Wang" cluster definition. For general and 
historical discussions of the Ising droplet definition, see Refs. ||146| , |147|| , and for discus- 
sions in the context of metastability, see Refs. [104-106, 145]. 

A serious problem of a numerical nature is the difficulty of accurately locating the 
coexistence curve in a model in which its position is not given by symmetry, as it is for 
the Ising or binary lattice-gas model. (Even in more complicated lattice-gas models this 

The extremely strong dependence of the nucleation rate 



is a nontrivial problem [ |14£ 
on the distance from the coexistence curve will cause even a small error in its location 
to produce large errors in numerical estimates of the droplet free energy and the pre- 
exponential factor. 

Recent large-scale computer simulations are discussed in Refs. 



30 , and further 



reviews of progress in the theory of metastability in fluids can be found in Refs. [22-24, 
31]. 



5 Ising Models with Weak Long-Range Forces 

In the limit of infinite interaction range, ferromagnetic Ising models with weak, long-range 
interactions (WLRF models) have been rigorously shown to possess infinitely long-lived 
metastable phases |58], [59|, |149|| . These phases lie on a hysteresis (or van der Waals) loop 
that reaches beyond the Maxwell construction in the phase diagram, and whose critical 
points correspond to sharp spinodals beyond which the metastable phases disappear. In 
this respect a WLRF model is similar to a mean-field theory || |||, p9|1 , except that in a 
mean-field theory the only allowed form of fluctuation is a spatially uniform change of the 
order parameter. The free-energy cost of such a change is extensive in the system volume 
and the mean-field approximation therefore predicts that metastable lifetimes 



135 



are infinite in the thermodynamic limit. 

Since the behavior of a number of physical systems, including superconductors and 
long-chain polymer mixtures, are often well described by WLRF models, there has been 
an interest in studying metastability in these models. Analytic approaches include the 
early field-theory treatment of a model fluid with conserved order parameter by Cahn and 
Hilliard 1T5L IT 



, the Fokker-Planck equation |89| , |135| , |150| , renormalization-group analyses 
p.51| , |152|| , and arguments from random long-range bond percolation ||153| , |154| , in addition 



to analytic continuation of the free-energy density [83, 92-95]. Numerical approaches 



10 



include MC simulations [129, 131, 134, 135, 153-155], traditional TM calculations [[15 
and calculations of |Im/| by the CTM method described in Sec. P^| [p4] ^5] [13§, p3fl . 



As pointed out by Klein and coworkers [92, 93, 153-155], the critical droplet in a 
WLRF model near the spinodal is a ramified structure whose local order parameter is 
only slightly different from that of the metastable phase. (Note the contrast with the SRF 
case, in which the critical droplet is compact. See Sec. |BTT| .) Only during the supercritical 
growth phase does the droplet center compactify, leading to a measurable change in the 
global order parameter. As a consequence, it is difficult to determine the nucleation 
time accurately from MC simulations by monitoring the order parameter or, equivalently, 
the nonequilibrium relaxation function. Therefore, several new techniques have been 
developed to analyze MC data for metastable decay in WLRF systems, including the 
analysis of recrossing events ||129| , monitoring the number of spins in the largest cluster 
[ |155|] , and intervention techniques [|154j] . Some of the recent work in this area has been 



reviewed in Ref. 11145 



The difficulties associated with accurately determining the nucleation rate in WLRF 
models make them particularly attractive candidates for study by alternative numerical 
techniques, such as the CTM method. A simple WLRF system suitable for CTM inves- 
tigation is the Quasi-One-Dimensional Ising (Q1D1) model, introduced by Novotny et al. 
15(ii|| . In its simplest form it is a chain of L layers, each of which contains N Ising spins, 



si tn =±l. Each spin interacts ferromagnetically with each of the 2N spins in the adjacent 
layers with coupling J/N>0 (these are the nonzero Jjy in Eq. ([!]) for this model), and 
with an external field H. In terms of the single-layer magnetizations, mi=N~ 1 J2n=i s i,n, 
the Hamiltonian is 



H = -NY\Jmimi+i + Hmi) . (13) 

i=i 



5.1 Nucleation 



In a WLRF system like the Q1DI model, where droplet "interfaces" are difficult to define, 
the onset of nucleation is detected by a rapid increase in the magnetization of some layer 
to within a small neighborhood of the equilibrium value. The physical picture is that 
of a bell-shaped magnetization profile that quickly develops into a radially expanding 
front separating the two competing phases. It was shown by Gorman et al. |94| that the 



free-energy cost of nucleation for the Q1DI model can be approximated by mapping the 
free-energy density functional of the model to a Ginzburg-Landau form and solving the 
resulting Euler-Lagrange equation. A 3 expansion of the potential around the exactly 



known |33|, |135|| mean- field spinodal field H s , in which X = \H\ — \H S \ —>■ , gives the 



following free-energy cost for a <i-dimensional WLRF system with force range R p2| , |9~3]j : 



F c = A(T)R d \X\ (6 ' d)/4 [l + 0(X 1/2 )] , (14) 

where A(T) is a nonuniversal function of T. For the Q1DI model R=N and d=l [ 152 , |156 



By solving the Euler-Lagrange equation numerically, one can also take into account the 
corrections to the expansion, as well as discrete lattice effects [J94fl . 

The prefactor to the Van't Hoff-Arrhenius term in |Im/| can be calculated by ex- 
panding the free-energy density functional about the stationary points corresponding to 
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the metastable phase and the critical droplet. The determinants of the two resulting 



Schrodinger operators give [M 



\lmf\ = B(T){V'/V)R d/2 \\\ d{1 - d/8) [l + 0(A 1/2 )] e"^ , (15) 

where B(T) is a nonuniversal function of T, V is the system volume, and V is the volume 
of the subspace in which the droplet itself is free to move without a cost in free energy. 
For the Q1DI model, V /V=L/(NL)=N~ 1 0. 

From Eqs. (|T^D and ([T5| ) we see that for large N, unless H is extremely close to H s , 
the free-energy cost F c of surmounting the nucleation barrier is large, so the exponential 
factor sets the scale for the metastable lifetime. However, for small N, or for H^H S , the 
lifetime is more strongly dependent on the particulars of the dynamic and on the detailed 
structure of the saddle point in the free-energy functional. The finite-range-scaling of 
|Im/| near the spinodal is found by recasting Eqs. (|HD and (|15|) in terms of a scaling 
variable (=R M ^ e ' d ^\X\, giving 

|Im/| = (V'/V)R- W2)[d+3{d - 2)/{6 - d)] $(T,C) (16) 

with a scaling function $ dependent only on T and (. The dynamic prefactor k, taken 
from the lowest eigenvalue of the Schrodinger operator corresponding to the saddle point, 
scales as (A) 1 / 2 , but is independent of R. The finite-range scaling of the nucleation rate T 
is thus 

r = (^v'/V)R~ {d/2)[d+{3d ~ 2)/{6 ~ d)] g(T,Q (17) 

with a scaling function Q dependent only on T and £. 
5.2 Numerical Transfer-Matrix Results 



The CTM method outlined in Sec. |3.2| was applied by Gorman et al. [Q to Q1DI sys- 
tems of infinite length and finite cross section N, with N up to between 100 and 500. 
Figure [3] shows typical spectra of Re/ a and |Im/ a | plotted against if at a fixed tempera- 
ture. These spectra are compared with the analytically continued free energy in the limit 
iV— >oo, which corresponds to a Curie- Weiss mean-field result (thick, dashed curves). In 
each calculation, as demonstrated in the figure, a unique a was found for which Re/ a 
computed from Eq. ([12]) closely approximated the real part of the mean-field metastable 
free energy. Invariably, the same a produced the smallest nonzero value for |Im/ a |. As 
H was changed away from H s towards H=0, the value of the metastable |Im/ a | was 
exponentially suppressed. 

If the metastable f a represents the analytically continued free-energy density, then 
one expects In llm/^! to be a measure of the height of the nucleation barrier, as seen from 
Eq. fll5D . Extrapolated CTM estimates of the barrier height F c are compared in Fig. [|with 
the free-energy cost of the saddle-point solution to the Euler-Lagrange equation, which 
was obtained by numerical integration. The agreement is remarkable and extends over 
an extremely wide range of fields and temperatures, for most of which the lifetimes are 
too long for standard Monte Carlo techniques to be useful. Similar results were found in 
a study of a three-state WLRF model |)]| , even where two competing metastable states 



were present, suggesting that Langer's formula has a wider applicability than previous 
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studies have claimed J85J . These CTM results, as well as the MC results in Ref. [ |95|1 , 
are consistent with Ostwald's empirical "Law of Stages" (Gesetz der Umwandlungsstufen) 
[|19| . |157|| , whereby a metastable phase may decay via intermediate metastable states before 
reaching equilibrium. 



6 Ising Models with Short-Range Forces 

In contrast to the situation for WLRF models, metastable phases in systems with short- 
range forces (SRF models) eventually decay, even though their lifetimes may be many 
orders of magnitude larger than other characteristic timescales of the system [17]. Thus, 



whereas mean-field theory provides a qualitatively acceptable description of metastabil- 
ity for WLRF systems in the long-range limit, it gives a quite misleading picture for 
SRF systems. The following subsections are devoted to a discussion of metastability and 
metastable decay appropriate for SRF systems. 

As a prototype for the metastable dynamics of SRF systems, we consider in detail the 
decay of the magnetization in an impurity-free kinetic nearest-neighbor Ising ferromagnet 
in an unfavorable applied field. The critical point of this model is in the same static 



universality class as the liquid/vapor phase transition | 158 |. (However, the liquid/vapor 
transition is in a different dynamic universality class: that of Model H ||61|| .) The Hamil- 
tonian is obtained from Eq. ([J) by setting Jij=J if i and j are nearest-neighbor sites and 
Jjj=0 otherwise, explicitly yielding 



H 



■ J Y, s i s j- H Y,' 



(id) 



where J2u,j) an d X^ run over ai l nearest-neighbor pairs and over all sites on a <i-dimensional 
hypercubic lattice of volume L d , respectively |159 [. 

In a typical MC study of metastable decay in this model one starts from the metastable 
phase and follows the relaxation of the magnetization, which is closely related to Binder's 



droplet size distribution | 

is 4> s (t) = (m ms -m(t))/(m ms -m s ) 



nonequilibrium relaxation functions ||119| and is directly obtained as an average over the 

123| , |132|| . The volume fraction of stable phase at time t 
where m s and m ms are the bulk equilibrium and 
metastable magnetizations, respectively. The metastable lifetime is typically estimated 
as the mean-first-passage-time for 4> s (t) to a preset value. Monte Carlo studies using this 
or similar methodology have been performed in two [37, 104, 108, 110, 111, 120-129], three 
[38, 105, 130-132], and higher |133|| dimensions. Several of these studies were analysed in 
terms of droplet theory, establishing general agreement between theory and simulations. 
A potentially more accurate, but also more computationally intensive, method to estimate 
the lifetimes could use the recrossing-event distribution discussed by Paul and Heermann 
[|129|| to determine a field-dependent cutoff value for <p s , as suggested by Rikvold et al. 

ra. 

Metastability in the two-dimensional nearest-neighbor Ising ferromagnet has also been 



studied by traditional TM methods by McCraw, Schulman, and Privman ||127j [I i 
and by the CTM method by Giinther, Rikvold, and Novotny |T3S|, [T3U], |g§ . 



141 
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6.1 Nucleation 

To obtain quantitative comparisons between numerical results and the droplet-based nu- 
cleation theory, one must calculate explicitly the field-theoretical expression for the nu- 
cleation rate, Eq. ([?]). This requires as input the free energy of a compact critical droplet 
F C (T, H), as obtained by CNT ff|, |90|. (Note the contrast with the WLRF case, in which 
the critical droplet is ramified. See Sec. [|.) Sufficiently far below the critical tempera- 



ture we can calculate this by a standard droplet-theory argument (see, e.g., [18|, [72|, |73|| ) 
that is modified to consider the nonspherical droplets that appear at low T because of 
the anisotropy of the surface tension [108, 138, 139, 144, 160-163]. The free energy of 
a (i-dimensional droplet of radius R (defined as half the extent of the droplet along a 
primitive lattice vector) and volume QdR d is 

F(R) = nf- 1)fd B^- x t - \H\AmQ d R d . (19) 

The quantity S is a temperature- and, in principle, field-dependent proportionality factor 
which relates the surface contribution to F(R) with the droplet volume ||162| , |164j| . The 



difference in bulk free-energy density between the metastable and stable states is \H\Am, 
which takes some account of droplet nesting through the magnetization difference Am 



60] , |165| . Maximizing F(R) yields the critical radius, 



flc(T,ff) = ^y^ (20) 

\H Am 



where ao = fl/(dCl d ) is the surface tension along a primitive lattice vector [162], and the 
free-energy cost of a critical droplet, 




F ^v = [m^) UJ ■ (21) 

In addition to R c , the critical droplet is characterized by other degrees of freedom, in- 
cluding the critical growth mode, droplet translations, and deformations represented by 
capillary waves on the droplet surface. The field-theoretical expression for the nucleation 
rate, Eq. ([?[), properly accounts for the effects of these additional degrees of freedom, and 
it is obtained explicitly by a saddle-point calculation that yields [76-78, 90] 

T(T,H) = A iT)\ H \b+o e -PF c {T,H){i + om) = A(T) |#|Hc e -(/E/|H|^)(i+0(fl*)) (22) 

with 



'd-iV- 1 




(23) 



y Am J 
The quantity A(T) is a nonuniversal function of the temperature only, 

f (3-d)d/2 forl<d<5,d^3 ( . 

b ~\ -7/3 for d=3 {2A) 

is a universal exponent related to excitations on the surface of the critical droplet |7(| [H| . 
and c gives the H dependence of the kinetic prefactor n |77|, ff8| . The kinetic prefactor is 
the only part of T(T, H) that may depend explicitly on the specific dynamic. 
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The Becker-Doring cluster dynamics is defined in terms of a master equation for the 
probability distribution q of /-particle clusters. The original theory only allows I to 
change by ±1 [3 1], but later modifications also allow cluster coagulation and fragmentation 
[|122|| . The CNT result for the exponential part of T(T,H) can be obtained from this 
approach. The existence and uniqueness of a solution q(£) that is metastable in the sense 
of the Penrose-Lebowitz-Sewell criteria have been proven both for the original dynamic 
[|166| , |167|1 and for the coagulation-fragmentation generalization [|168| . However, describing 



the clusters only in terms of their particle numbers does not suffice to obtain the pre- 
exponential factor in Eq. (0) for a particular system. This was demonstrated explicitly by 
Binder and Stauffer ||125|| , who obtained a result formally analogous to Eq. (|22"D by using 
a droplet model in which the individual clusters were characterized by several coordinates 
in addition to the particle number I. 

For d=2, there is substantial numerical evidence that b=l, as predicted by Eq. (^4j). 
This is obtained from calculations that do not involve the dynamics, such as analyses of 
series expansions ||160| , |169| , |1 70|| and transfer-matrix calculations ||138| , |139| , |144j| . These 
studies, as well as MC work [108, 171-173], also indicate that the free-energy cost of the 
critical droplet is given by Eq. fl2"T|) with the zero-field equilibrium values for £ and Am. 
We therefore adopt the notations S=S(T), Am=2m s (T) , and S=S(T) to emphasize the 
lack of field dependence in these quantities. The quantity S(T) can be obtained with ar- 
bitrary numerical precision by combining a Wulff construction with the exact, anisotropic 
zero- field surface tension [|162|| , and m s (T) is obtained from the exact Onsager-Yang equa- 
tion ||1 74| . These general results, that the surface free energy and bulk magnetization of 



compact critical droplets are determined by the zero-field equilibrium surface tension and 
magnetization, respectively, are also supported by MC studies of nucleation rates in three 
dimensions |D§ [ISO], P2] , 

For dynamics that can be described by a Fokker-Planck equation, it is expected that 
the kinetic prefactor is proportional to R~ 2 |77|, [78], [Kj, which by Eq. ( p0[) yields c=2. 



This value has been confirmed numerically for the Metropolis algorithm with updates at 
randomly chosen sites, but it does not appear to apply if the sites are chosen sequentially 

The numerical results cited above indicate that CNT, modified by the post-classical 
results for the prefactor exponents, gives very good agreement with the observed behavior 
of kinetic Ising and lattice-gas models, even moderately far away from the coexistence 
curve. A rough estimate for the field strength beyond which the simple droplet theory 
discussed here should break down (or at least become suspect) is obtained by requiring 
that 2R C (H,T)>1 ||139| . The resulting crossover field is sometimes called the "mean- field 

HTD8L 



spinodal point", or MFSP 



IT01 



but it is not identical to the sharp spinodal found 
when the Ising ferromagnet is treated in the mean-field approximation. The MFSP is 
located at Hmfsp{T) = 2(d — l)<7o(T)/2m s (T). The field region beyond this limit we call 
the strong-field region. It will not be discussed further here, but we hope to return to it 
in the future. 



6.2 Growth and the KJMA Theory 

The reasoning behind the Kolmogorov-Johnson-Mehl-Avrami (KJMA) theory is quite 
simple [96-98, 175-178]. Critical droplets are assumed to nucleate in the metastable 
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phase and subsequently grow without substantial deformation ||179|1 . (Recent discussions 
relevant to the limitations of the latter assumption can be found in Refs. [180-183].) The 
nucleation rate T may either be constant (homogeneous nucleation), or all the nuclei may 
already be present at t=0 (heterogeneous nucleation |33)). The KJMA theory is simple 
to work out for both cases |33], |68| , but only the former will be pursued here. 

The radial growth velocity, which approaches a constant limit for large droplets 
obtained in an "Allen-Cahn" approximation [72, 73, 175-178, 184-186] as 



is 



v± = (d-l)u (R^-R- 1 ) ^ (d-l)uR- 1 = v , 



(25) 



where the coefficient v depends on the details of the kinetics. (To avoid confusion, we 
note that the usual growth law for spinodal decomposition in systems with nonconserved 
order parameter, i?~t 1//2 , results from setting i?~ 1 =0 (i.e., H=0) in Eq. Q2"5p. The i?~t 
growth law obtained in the case of metastable decay is a consequence of the difference 
in bulk free energy between the two phases, which acts as a driving force for the growth 



184, 185 



If we set v±=Vo, neglect the volumes of the critical droplets, and consider an uncor- 
related "ideal gas" of freely overlapping domains of the stable phase, then the volume 
fraction of stable phase at time t is 



Mt) 



exp 



-m d v$ f\t-s) d ds 

Jo 



exp 



a, 



d+1 u 



t\ d+1 ' 



(26) 



where the integration variable s is the time at which a particular supercritical droplet was 
nucleated. This relation is known as Avrami's law. The timescale to is the characteristic 
time for collisionless growth and sets the basic timescale for the decay of the metastable 
phase. By performing the integration in Eq. (p6|), one sees that to is given by 



d-p\ — 



t (T,#) = Kn 



d+1 



B(T)\H\ 



b+c + d 

d + 1 exp 



1 /3H(T) 
d + llHl^ 1 



(27) 



where B(T) is a nonuniversal function of T. The second equality in Eq. (|27| ) was obtained 
by using Eq. (|22]) for the homogeneous nucleation rate T(T,H), neglecting for simplicity 
the correction term in the exponential. By comparing Eq. (^) for t with Eq. ( p2| ) for T, 
we notice that, apart from the pre-exponential factors, the characteristic time t is simply 
given by the nucleation rate to the power — l/(d+l) [|132| . 

Associated with t is the characteristic lengthscale for collisionless growth, which we 
loosely call the mean droplet separation. It is given by 



R (T,H) = v t = C(T)\H\ 



b+c-l 

d + 1 exp 



1 /?S(T) 
d+l\H\ d - 1 



(2c< 



where C(T) is a nonuniversal function of T. The lengthscale Rq can be very large for 



\H\ 



•oo as 



weak fields, and we note that, even though the critical droplet radius R c 
\H\ — >0, R c /Ro^O in the same limit. 

The KJMA theory gives a good approximation for sufficiently small <ft s (well below 
the percolation limit |147j| ) that droplet correlations do not significantly alter the growth, 
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and it quickly became popular for analyzing experimental results on metastable decay in 
a number of areas. These include situations in which d is smaller than the dimension of 
the physical space and must be interpreted as the dimension of the subspace in which 



the droplets grow p8| . However, the theory does not seem to have attracted sustained 
attention from theorists until the 1980 's, when Sekimoto derived exact expressions for the 
space-time correlation function and the structure factor for the KJMA process [175-178]. 
These results were later generalized to infinitely [|187|1 and multiply |188|1 degenerate equi- 
librium phases, and mean-field results for the multiply degenerate case with multiple nu- 
cleation and growth rates have been obtained, both with homogeneous and heterogeneous 
nucleation ||189| , 190|| . Applications of Eq. (|26|) have recently been made to theoretical 
studies of metastable decay in long-range interaction models and ferroelectrics [35-38], to 
nanometer-sized ferromagnetic particles |49] , and to kinetic Ising models |37|, |38], [E|, |108 



Several of these papers contain kinetic MC simulations [37, [H|, ^, |108| , |187| |189j , |190 



6.3 Finite-Size Effects 

In an infinitely large system, the number of critical and supercritical droplets is of course 
infinite, and the decay of the order parameter is well described by Avrami's law, Eq. (ffijj). 
In this section we consider how finite system size modifies the KJMA picture. Although 
finite-size scaling has been extremely useful in the study of equilibrium critical phenomena 
(see, e.g., Ref. ||191|| ), systematic finite-size analysis of metastable decay seems to have 
been performed only recently [[36], |37|, |108j , |110|| . The discussion below follows closely that 
of Ref. [|108j| . For simplicity we retain our exclusive focus on hypercubic systems of volume 
L d with periodic boundary conditions. 

For temperatures well below the critical temperature T c and fields inside the MFSP, 
the correlation lengths in both the stable and the metastable phase are microscopic. We 
are then left to consider the interplay between three lengths: the system size L, the mean 
droplet separation i? , and the critical radius R c . 

In the large- L limit, where 



L > R > R c 



(29) 



the system can be approximately partitioned into (L/Ro) d ^> 1 cells of volume R d . Each 
cell decays in an independent Poisson process of rate R d T = tQ 1 . The volume fraction is 
then self-averaging, so that Eq. (|26|) can be inverted to yield the average time it takes for 
the volume fraction of the equilibrium phase to increase to a given value <f) s , 



(t(<f> B ))*to(T,H) 



The relative standard deviation of tid*) is 



d+1 



ln(l- 



i 

d+l 



(30) 






(Ro/L)*p, 



(31) 



where p~l is the relative standard deviation of a single Poisson process (not to be confused 
with the lattice-gas density defined in Eq. (||)). 
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The regime characterized by r ^C 1 has been termed "the deterministic region" ||108 , 

110 |. It is subdivided into the strong-field region mentioned in Sec. |6J], and a region 

characterized by a finite density of growing droplets, which we call the multi- droplet region 

108|| . Observations of the deterministic region in MC simulations are also indicated in 

!^, m |105| , |108| , |122| , |126| , |132|| . The characteristic absence of L-dependence in 



Refs. 

the multi-droplet regime was noted by Binder and Muller-Krumbhaar, who also derived 



equations equivalent to Eqs. 



and (M 



|122| . Although the main emphasis was on the 
nucleation process, the multi-droplet picture was also implied in Langer's work [76-78]. 
For smaller L, so that 

# > L > R c , (32) 

the random nucleation of a single critical droplet in a Poisson process of rate L d T is 
the rate-determining step. This is followed by relatively rapid growth, until this droplet 
occupies the entire system after an additional time much shorter than the average waiting 
time before a second droplet nucleates. Therefore, the characteristic lifetime becomes 



<*(&)>« (L d T(T,H) 



L- d [A(T)]^\H\ 



-(b+c) 



exp 



/3~(T) 



\H 



d-l 



(33) 



In this case r ~ 1, and (t(<p B )) depends only weakly on the threshold <fi a . This single-droplet 
region is part of "the stochastic region" identified in Ref. ||1 10[ , and it was detected in 



MC simulations in Refs. 



The crossover between the deterministic and stochastic regimes is determined by the 
condition LccRq with a proportionality constant of order unity. We identify the crossover 
field with the "dynamic spinodal point" (DSP) introduced in Ref. [f LlOj , and in the limit 
H—>0 we explicitly obtain from Eq. (|28|) 



tfi 



DSP 



(InL)- 



i 

' d-\ 



(34) 



This crossover field was observed in MC simulations reported in Refs. [37], 0, £|9], |105| , 



108| , |126| , |132| . Following Refs. [|108| , |110|1 , we estimate i^osp as the field where the rel- 
ative standard deviation for the lifetime is r=l/2. We emphasize that, although -/^dsp 
vanishes as L — > oo, the approach to zero is exceedingly slow, especially for d=3 and 
above. Therefore, i^DSP may well be measurably different from zero for systems that are 
definitely macroscopic as far as their equilibrium properties are concerned. (As an illus- 
tration, increasing L from 100 to 10 10 for d=3 decreases the leading term in i^osp only 
to approximately one-half of its original value!) 
Finally, we consider the small- L limit, 



i?o > #c > L 



(35) 



In this case the volume term can be neglected in Eq. (|T9|), and the free-energy cost of 
a droplet occupying a volume fraction </> s = V(R)/L d is -F(0 S ) ck L d ~ l (f)^ d ~ 1 ^' d Yi(T) with 
a proportionality constant between 1 and 1/d, so that the first-passage time to a given 
4> s is independent of H and diverges exponentially with L d ~ x . Since the dynamics in 
this region of extremely weak fields or extremely small systems is similar to that on the 
coexistence line, H=0 | 192 |, we call it the coexistence region. The crossover field between 
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the coexistence and single-droplet regions, called "the thermodynamic spinodal point" 
(THSP) in Ref. [ITO], is determined for a given S by Qd{Rc/L) d = (p s , which yields 



_ 1 (d-l)S(T) 
FTHSP " L^ 2dm s (T) ■ (36) 



This crossover field was observed in MC simulations reported in Refs. [^, |108| , |127 [ 



In summary, by comparing the characteristic lengths Rq and R c with the lattice con- 
stant and the system size L, one can identify four different field regions, in which the decay 
proceeds through different excitations. In order of increasingly strong unfavorable field 
\H\, these are the "coexistence region," characterized by subcritical fluctuations on the 
scale of the system volume; the "single-droplet region," characterized by decay via a single 
critical droplet; the "mult i- droplet region," characterized by decay via a finite density of 
droplets; and the "strong-field region," in which the droplet picture is inappropriate. The 
crossover fields between these regions, 



-l 1 



[Htksp ~ £ ] < [#dsp ~ (In L) a=r] < [#mfsp ~ L v ] , (37) 

are accurately predicted by droplet theory. The different regions and crossover fields are 
illustrated in Figs. [3)-||. 

In Fig. Q are shown average metastable lifetimes for LxL square Ising ferromagnets 
with L=128 and 720 at T=0.8T C . The data points were obtained by the Metropolis 



algorithm with updates at randomly chosen sites |108|1 and have here been extrapolated 
past i^THSP f° r £=128 (given by Eq. Q3T5D), using Eq. Q5BD for the lifetime in the single- 
droplet region. The four field regions can clearly be distinguished. 

An alternative view of the information contained in Fig. ^ is found in Fig. |J which 
shows as a function of L the field H sw at which (t(0 s =l/2)) = r, plotted for two different 



values of r at T=0.8T C j49[. This figure corresponds to a contour plot of data like those 
shown in Fig. |3|. In accordance with usage in the experimental literature on small magnetic 
particles [41-48], we call H sw the "switching field." For qualitative comparison we have 



also included in Fig. |4] data digitized from Fig. 5 of Ref. [45|| , which shows the effective 



switching field (corrected for the demagnetization field) vs. the particle diameter for single- 
domain ferromagnetic barium ferrite particles, measured by magnetic-force microscopy 
at room temperature. Considering the different dimensionalities of the model and the 
experimental system, and that no particular effort was made to fit the parameters in 
the Ising model to the experiments, we find the similarity between the simulated and 
the experimental switching fields striking. This qualitative agreement may indicate that 
the decay of the metastable magnetization state in the barium ferrite particles proceeds 
through similar nucleation and growth mechanisms as in the Ising model, and it should 
be relevant to the current debate over the magnetization reversal mode in single-domain 
ferromagnetic particles [42-48]. 

In both the multi-droplet and the single-droplet regions, the metastable lifetime (de- 
termined by Eqs. (|27| ) and (^), respectively) has the form of an exponential in l/\H\ d ~ l 
multiplied by a power-law prefactor in \H\. In both regions the derivative of ln(t(0 s )) 
with respect to l/\H\ d ~ 1 can therefore be written as 

_ dln(t(0 s )) xio-id-i, A /ooa 

cff = d(l/|ff| rf - 1 ) ' ' ' ( ' 
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with A = (b+c)/(d— 1) and A = /3H(T) in the single-droplet region, and X = (b+c+d)/(d 2 — 1) 
and A=/?E(T)/(d+l) in the mult i- droplet region. The quantity A e g, calculated from the 
MC data in Fig. |3], is shown in Fig. |5|. The dashed straight lines correspond to Eq. fl38|) 
with the theoretically expected exponent b+c=3 and the numerically exact S(0.8T C ) ||162 
(see discussion in Sec. |6.1|). The data points for both system sizes follow the lower of the 



two lines in the weak-field part of the multi-droplet region. The steep rise in A e g expected 
near -£/dsp is seen for both systems. However, only the smaller one penetrates into the 
single-droplet region in the field range for which data could be obtained with a reasonable 
amount of computer time. A detailed statistical analysis is given in Ref. ||108| . 

Due to the dramatic increase in the metastable lifetime as T is lowered, numerical 
data confirming the theoretical predictions at lower temperatures must be obtained by 
other techniques, such as the CTM method or one of the new MC methods discussed in 
Sec. |||. In Fig. |6] we show the temperature dependence of i^osp f° r LxL systems with 
L=24 and 240, as obtained by the new MC with absorbing Markov Chains (MCAMC) 
method [113-115], together with our analytic estimate for H MFSP . The slow decay of 
i^DSP with L predicted by Eq. (|34]) can be seen. Also note the dramatic widening of the 
single-droplet region as T is lowered for a system of fixed size, in agreement with recent 
exact predictions |106| , |1 93 . 



194. 



Confirmation that the quantity 5(T) is given by its zero- field equilibrium value is 
given in Fig. [7] (after Ref. ||195| ), which shows estimates of S(T) for T between 0.17T C 
(T/J=0.4) and 0.8T C based on the CTM method |3|, p9| , MCAMC simulations [|TT5j 
and standard Metropolis MC simulations 



108 



The CTM estimates were obtained by 
fitting the logarithmic derivative of the metastable |Im/ Q | to Eq. (^) with 6=1, c=0, and 
the correction term from Eq. (f22|) included. The numerical estimates for H(T) obtained by 
the different methods agree very well with each other, as well as with the exact equilibrium 
result. 



6.4 Discrete-Lattice Effects 

Recently, a series of rigorous papers [106, 193, 194, 196-198] have appeared that dis- 
cuss effects of the lattice discreteness on the metastable lifetimes in the limit T—>0. The 
approach is to consider a birth-death process describing the time dependence of the fluc- 
tuations corresponding to subcritical droplets in the two-dimensional Ising ferromagnet 
with single-spin-flip Metropolis |193| , |194j , |196| , |197|| or modified Swendsen-Wang ||106 



dynamics. For isotropic interactions | ]106| , |193| , |194| , |196| ] and \H\/J<2 it is shown that 
the metastable phase almost certainly decays through a single "proto-critical" droplet of 
spins pointing parallel to the field. This droplet is shaped like a rectangle of Z c x(Z c — 1) 
overturned spins, with one additional overturned spin attached as a "knob" to one of 
its long sides. The length l c =\2J /\H\\ is the smallest integer larger than 2J/|i?|, where 
2 J J | H | = lim-r^o 2-R c (T, H) is the zero-temperature limit of the diameter of the critical 
droplet. In the zero-temperature limit the average metastable lifetime is found to be a 
piecewise linear function in \H\ which diverges asymptotically as 1/|-H"| for small \H\. It 
is given by 



lim£;BTln(L 2 (£(0 s = l))) = 8Jl c -(^-/ c +l) 2|if| 



T->0 



8J 2 /\H\+U-0(\H\) . 



(39) 
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These results are significant for two reasons. First, they provide an explicit dynamical 
justification for the existence of a critical droplet size in kinetic Ising models. Second, 
they generalize the standard continuum droplet theory outlined in Sec. |6.1| to consider the 
discreteness of the lattice. Using the exact zero-temperature value, S(T=0) = 8 J, one sees 
that Eq. (|39|) is consistent to leading order in J/\H\ with the continuum droplet-theory 
result for T, Eqs. © and (H). 



The temperatures at which these discrete-lattice results are expected to be valid are 
so low as to be inaccessible with standard MC techniques. However, they are readily ac- 
cessible, both with the MCAMC method [113-115] and with the CTM method |13|, |T39] . 
Results obtained with these methods are shown in Fig. || together with the derivative of 
the first line in Eq. (|39"D with respect to 1/\H\ and the corresponding quantity in contin- 
uum droplet theory, A; B TA eff from Eq. fl3~8l). The discrete-lattice results are represented by 
the series of solid, parabolic arcs, and the continuum results for the single-droplet region 
are represented by the two straight lines. The dashed line corresponds to b=l and c=2, 
which is appropriate for the MC results, whereas the dotted line corresponds to 6=1 and 
c=0 and is appropriate for the CTM results which do not contain a kinetic prefactor. 
For relatively strong fields, the MC data agree reasonably well with the discrete-lattice 
predictions. For weaker fields, where the critical droplets become larger, the oscillations 
caused by the lattice discreteness become less pronounced, and the MC data points ap- 
pear to approach the continuum result. The CTM method allows a closer approach to 
H=0 than the MC, but the deviations from the continuum result are larger than for the 
MC at the same field. A detailed comparison of the manner in which the results obtained 
by the two methods approach their respective continuum limits, represented by the two 
straight lines which intersect at the common zero-field limit £, has yet to be performed. It 
is likely to involve both the difference between the geometries of the two systems studied 
(square for MC and infinite strip for CTM) and the fact that the CTM quantities can 
be evaluated only at particular fields determined by the lobe structure of the metastable 
IIhi/q,! [ [L38| , 139 1, which is similar to the one shown for the WLRF Q1DI model in Fig. [I]. 



The theoretical and numerical results discussed in this subsection raise a number of 
questions that should be answerable using the new MC and TM algorithms. In particular 
it is important to understand how the lifetimes cross over to the well confirmed results of 
continuum nucleation theory at higher temperatures and weaker fields, and also how the 
lifetimes and the crossover to continuum theory are affected by a change from Metropolis 
to Glauber dynamics. Furthermore, some rather dramatic effects have been predicted for 



the case of anisotropic interactions as T^O |197|| . No sign of these effects were found in 



a recent study by the CTM method |144fl , and a further investigation by the MCAMC 
method is in progress. 

7 Summary and Discussion 

In this review we have presented some aspects of the current knowledge regarding the 
mechanisms and rates of decay of metastable phases in statistical-mechanical systems. 
With its record of 270 years of published research in a variety of basic and applied contexts 
following Fahrenheit's 1724 paper [||], this is indeed a venerable field of inquiry. Yet we 
believe it is a sign of continued vigor and relevance that nearly half of the references we 
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cite have appeared within the last decade. 

It is now well understood that metastable decay is a kinetic phenomenon, whose rate 
depends on the interplay between the (homogeneous and/or heterogeneous) nucleation of 
critical droplets of the equilibrium phase and on the subsequent growth of the supercritical 
droplets. Formally, the nucleation rate is given by the Van't Hoff-Arrhenius relation of 
classical nucleation theory, Eq. (H), or the field-theoretical relation, Eq. (|7]) [76-78], and 
the simultaneous nucleation and growth give rise to "Avrami's law", Eq. (^) [96-98]. 
However, both the free energy F c of the critical droplet and the pre-exponential factor V 
in Eq. (^|) (both of which are contained in the imaginary part of the metastable free energy, 
Im/, in Eq. (|7|)) depend crucially on the geometry of the critical droplet. Whereas the 
shape of the critical droplet determines the exponential factor in Eq. (|6]), the prefactor is 
determined by the fluctuations around the critical droplet shape. The latter is evaluated 
in the "post-classical" field-theoretical approach by a steepest-descent calculation around 
a saddle point representing the critical droplet. 

If the order parameter is constrained to be uniform over the entire system, so that 
F c becomes extensive in the system size, one recovers the mean-field picture of van der 
Waals @, Maxwell 0, and Neel [|39fl . In this picture the metastable phases are infinitely 
long-lived in the limit of infinite system size, and the limit of metastability is marked by 
a sharp spinodal line. 

A different class of systems are the weak-long-range-force ( WLRF) models discussed 
in Sec. |5|, which may be useful models for, e.g., superconductors, long-chain polymers, and 
systems with elastic interactions. These models also have infinitely long-lived metastabil- 



ity and a sharp spinodal in the limit of infinite interaction range [|5£| [5^, |149|| . However, 
the critical droplet is a ramified object [92, 93, 153-155] whose free energy is given in 
terms of the interaction range and the distance from the spinodal by Eq. ([TJP |9~2", [£fl , 



and the resulting nucleation rate is given by Eq. (|TT|). The agreement between these ana- 
lytical results and the numerical constrained-transfer-matrix (CTM) method is shown in 
Fig. |)4|], which illustrates the dramatic increase in metastable lifetime as the distance 



from the mean-field spinodal increases. 

For short-range-force (SRF) models, which are useful to represent, e.g., anisotropic 
magnets governed by exchange interactions, and which belong to the same static univer- 
sality class as the liquid/vapor phase transition, the situation is again different. In this 
case the critical droplet is a compact object with a radius R c given by Eq. (|20D , and F c 
is given by Eq. fl2~T|) . The nucleation rate from the field-theoretical saddle-point calcula- 
tion is given in Eqs. (|22|) -(p4|) [76-78, 90], and the effects of simultaneous nucleation and 
growth are reflected in the explicit form of Avrami's law, Eqs. ( ^6|) and ( ^7|) [96-98]. In 
contrast to mean-field and WLRF systems, SRF systems do not have a sharp spinodal. 

Finite-size effects represent an aspect of metastable kinetics which has attracted in- 
creased attention in recent years |36], |37|, |108| , |110|| . For SRF systems, these effects give rise 



to different dependences of the metastable lifetime on applied field (or chemical potential, 
pressure, or supersaturation) and system size for different values of these parameters. 
These regions are separated by size dependent crossover fields: the thermodynamic spin- 
odal field i^THSP given in Eq. (|36|), the dynamic spinodal field i^osp given in Eq. (|34|), 
and the size- independent "mean-field" spinodal field H MFSP . F° r -^thsp< 1-^1 <-^dsp the 
decay occurs through a single droplet, and the mean lifetime is given by the inverse nu- 
cleation rate through Eq. (|33|). In contrast, for H DSP < \H\ <H MFSP the decay involves a 
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nonzero density of simultaneously nucleating and growing droplets, and the mean lifetime 
is given by Avrami's law through Eqs. fl2"T|) and fl30|). The detailed agreement between 
these analytical results and Monte Carlo (MC) simulations is illustrated in Fig. |^, which 
shows the dramatic dependence on l/|iJ| d_1 of the logarithm of the average metastable 
lifetime for two-dimensional systems of different sizes. Similar agreement has also been 
demonstrated for three-dimensional systems [130-132], and a rigorous justification for 
the nucleation-theory results has been obtained in two dimensions at very low tempera- 
tures [106, 193, 194, 196-198]. Further agreement between analytic theory and simulation 
results is illustrated in Figs. |, 0, and |8|. 

A number of questions regarding metastable decay still remain to be answered. Prob- 
ably some of the most difficult ones are related to fluids, which we discussed briefly in 
Sec. §. In particular we would like to see developed a droplet definition as unambiguous 
as the "Fortuin-Kastelyn-Coniglio-Klein-Swendsen-Wang" definition [104-106, 145-147] 
currently used for Ising and lattice-gas systems. 

Some of the most exciting potential for both fundamental and technological progress 
may lie in combining modern "atomic engineering" techniques, such as atomic and mag- 
netic force microscopies, with computer simulation. Using these experimental techniques, 
one could both build and study well-characterized metastable systems. Employing state- 
of-the-art computer hardware and algorithms, one could simulate these same systems 
with tens to billions [ |109| , |132|| of particles. This approach should provide the opportunity 
to study effects of size, boundary conditions, and impurities by a combination of theory, 
simulation, and experiment in a carefully controlled manner. A step in this direction 
was recently taken by Richards et al. |49|], some of whose results are shown together with 



experimental data from Ref. |45|] in Fig. f|. 

In summary, we think the field of metastability is a healthy 270- year-old with a long, 
exciting life ahead of it. 
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Figure 1: The branches of Re/ Q (top) and |Im/ a | (bottom) computed from constrained 
transfer matrices for a 35xoo Q1DI system at T=0.5T C and plotted versus H in units of 
j3~ l =T c =2J. The equilibrium branch of |Im/ Q | is identically zero and has been removed for 
clarity. The thick dashed lines represent the mean-field values for the real and imaginary 
parts of the analytically continued free energy. The dotted vertical line indicates the 
mean-field spinodal field H s . After Ref. P|. 



Figure 2: Extrapolated finite- iV CTM estimates for F c /N in Eq. (|TJ]), plotted in units of 
/5 C T 1 =T C =2J against \\\/H s = (H s — H)/H s for various temperatures T. The lines indicate 
the free-energy cost of forming the saddle-point solution of the Euler-Lagrange equation 
at the same temperatures and were obtained by numerical integration. After Ref. JJ|] . 



Figure 3: Metastable lifetimes for LxL two-dimensional Ising ferromagnets with L=128 
(o) and 720 (solid diamonds), shown on a logarithmic scale vs. J/\H\. The data were 
obtained by the Metropolis algorithm with updates at randomly chosen sites for T=0.8T C . 
The lifetimes, which are given in units of MC steps per spin (MCSS), were estimated 
as the mean-first-passage times to m=0.7 (0 S ~O.13) ||108f| . The solid curves represent 
extrapolations to \H\<Htyl$p for the smaller system, using Eq. ( J33|) for the lifetime in the 
single-droplet region. H^ SP is given by Eq. (|36| ) and is indicated by the arrow marked 
THSP(128). The other arrows mark the mean- field spinodal (MFSP) for both systems 
and the dynamic spinodal points (DSP) for L=128 and 720 respectively. The field- regions 
seen in the figure are, from left to right, the strong-field region, the multidroplet region, 
the single-droplet region, and the coexistence region (the latter for L=128 only). The 
dashed lines correspond to the asymptotic slopes in the single- droplet and multidroplet 
regions, (JEL and /3H/3, respectively. Data courtesy of S. W. Sides. 



Figure 4: The switching field H sw (L,t)/J (left-hand vertical axis) shown vs. L (bottom 
horizontal axis) for waiting times r=100 MCSS (o) and 1000 MCSS (x) in an LxL 
Ising ferromagnet at T=0.8T C . The data were obtained by MCAMC simulations, and the 
error bars, everywhere much smaller than the symbol size, are not shown. This figure 
corresponds to a contour plot of lifetime data like those shown in Fig. ^. The dotted 
curve represents i^THSP an d the dashed curve represents .£/dsp- To the left of .Hthsp 
is the coexistence region (CE), between the two spinodals is the single-droplet region 
(SD), and to the right of i^DSP is the multidroplet region (MD). Data courtesy of H. L. 



Richards and S. W. Sides [ISfl. For qualitative comparison we have also included data 
digitized from Fig. 5 of Ref. [45| (O with error bars), showing the effective switching field 



(right-hand vertical axis) vs. the particle diameter (top horizontal axis) for single-domain 
ferromagnetic barium ferrite particles. 
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Figure 5: Here are shown vs. \H\fJ estimated values of the quantity k^Th^/ J 2 from 
Eq. (p8|), based on the data from Fig. |[ The upper and lower dashed straight lines 
correspond to the theoretical result in the single-droplet and in the multidroplet region, 
respectively. Both were calculated using the theoretically expected exponent 6+c=3 and 
the numerically exact 5(0. 8T C )/ J 2 ?»0.918 915. The pairs of opposing arrows with error-bar 
"feathers" mark i^osp for L=128 (right pair) and 720 (left pair). See detailed discussion 



in Sec. 6.3. After Ref. 108 



Figure 6: The dynamic spinodal field, Hnsp/J vs. k-gT/J for LxL systems with L=24 
(o with error bars) and 240 (x with error bars), obtained by the MCAMC algorithm. 
The dotted curves are merely guides to the eye. Data courtesy of M. A. Novotny 
Also shown is our analytic estimate for the L-independent H MFSP /J (solid line). 



Figure 7: The quantity E(T)/J 2 as obtained from the CTM method (o) [|13|, [139], 
MCAMC simulations (solid diamond) |115|| , and Metropolis simulations (solid square) 
[ |108|| , shown vs. k^T/J. The solid curve corresponds to the zero-field equilibrium value 
of S(T) |162||. The corresponding critical-droplet shape interpolates between a square at 



T=0 and a circle at T=T C , both shown for the whole temperature range as dashed curves. 



See detailed discussion in Sec. |6T3| . After Ref. | |195|| 



Figure 8: Here are shown vs. \H\/J estimated values of the quantity k^TA.^/ J 2 from 
Eq. (|38|), based on metastable lifetimes for a 24x24 Ising ferromagnet at T=0.4J as 
obtained by the MCAMC method [|113| , |114|| (o with error bars), and on the metastable 
ln|Im/ a | for a 9xoo system at the same temperature, obtained by the CTM method ||139 1 
(O). The theoretically expected results from continuum droplet theory as given by the 



rhs of Eq. (|38|) are represented by straight lines. The dashed line corresponds to 6=1 and 
c=2, appropriate for the MC results, and the dotted line to 6=1 and c=0, appropriate 
for the CTM results, which do not contain a kinetic prefactor. Both lines intersect the 
y-axis at the numerically exact H(T=0.4J)/J 2 . An estimate of the effects of the lattice 
discreteness is given by the solid parabolic arcs, which represent the 1/|-H"| derivatives of 
the first line in Eq. (^). Data courtesy of M. A. Novotny and C. C. A. Giinther. After 
Ref. |IT4| . 
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